CAT plot for simulation

The population (dataset of interest) is 24 mice, 12 of each sex and 12 of each POE, with the same genetic composition as a result of crossing. Genes where at least three samples do not have 10 count were removed. Genes without at least one count for both alleles were removed. Genes with a significant sex, parent or interaction effect was removed.

Our sample size is 4 vs. 4. We select 8 subjects from our dataset, 4 for each sex and 4 for each POE. For the (i,j)th gene-sample in our real (sub)dataset, we simulate \(\text{LOGIT}_{ij} = \log\left(\frac{p}{1-p}\right)_{ij} = \alpha_i+\beta_i x\), where \(\alpha_i\) is the mean proportion of reads mapped to allele 1 for the ith gene, \(\beta_i \sim N(0,1)\) and \(x = (0,1,0,...,1)\). Then \(p_{ij} = \frac{1}{1+\exp(-\text{LOGIT}_{ij})}\). We then simulate \(Y_{ij}\) from \(\text{Betabin}(p_{ij}, n_{ij}, \phi_i)\) where \(n_{ij}\) is the total reads mapping to both alleles for gene-sample (i,j) and \(\phi_i\) is the overdispersion MLE of the ith gene. These will be our simulated allele counts.

Ash has a higher concordance than apeglm. Analysis last week also showed that ash has more extreme and frequent shrinkage, and has more accurate predictions on average as evaluated by mean absolute error.

Filtering

We wish to determine the filtering rule such that the MLE is most accurate with regard to a CAT plot. We look at three rules: removing genes where less than half the samples have a certain count (which we call the threshold), removing genes where less than all the samples have a certain count (threshold), and removing genes where the sum of counts across samples is less than a certain amount (threshold). For each rule, we look at various different thresholds, and calculate the “average MLE concordance” that results when filtering with that threshold. The average MLE concordance is just taking the concordance when looking at top 10 genes, top 20 genes, etc., and averaging. We chose the rule that has the best concordance, and if multiple rules have similar concordance, we chose the rule that has the smallest number of genes removed. Across all rules and thresholds, filtering out genes with total counts less than 610 is the best. Using this rule, MLE with filtering is neck-and-neck with ash in the CAT plot, and better than apeglm.

## Mean absolute error, no filter

## $MAEmle
## [1] 0.1556092
## 
## $MAEmap
## [1] 0.1446702
## 
## $MAEash
## [1] 0.143204

## MAE for apeglm-shrunk genes (shrinkage>0.1), no filter

## $MAEmle
## [1] 0.5114665
## 
## $MAEmap
## [1] 0.4195591
## 
## $MAEash
## [1] 0.405118

## MAE for ash-shrunk genes (shrinkage>0.1), no filter

## $MAEmle
## [1] 0.4763569
## 
## $MAEmap
## [1] 0.3850236
## 
## $MAEash
## [1] 0.3706626

## MAE for MLE with filter

## [1] 0.1012654

CAT plot for alternate simulation

We try simulating \(\beta_i\) from \(t_3/10\) so that we have mostly close-to-zero effects (but with some relatively large effects occasionally appearing). Furthermore, we simulate \(\phi_i\) from \(\text{Exp}(1/87)\) with probability \(0.5\) and 500 with probability \(0.5\). The distribution of overdispersion from the real data was approximated by a similar mixture distribution: one component was \(\text{Exp}(1/196)\) with 30% proportion, and the other component was a point mass at 500 with 70% proportion. Therefore, our simulated overdispersion will lead to nosier data.

Filtering for alternate simulation

In this case, no amount of filtering leads to concordance at top genes as good as apeglm and ash. Apeglm had an average concordance of close to 0.5.

5 vs. 5

We tried simulating standard normal effects in 5 vs. 5 sampling. Output is hidden as results are the same as the 4 vs. 4 normal effects case: ash gives more shrinkage, better estimates on average, and higher concordance. Filtering did not lead to universally better/worse concordance than ash (though it beat apeglm).

3 vs. 3

We simulated standard normal effects in a 3 vs. 3 design as well. Most output is hidden as results are similar to the 4 vs. 4 normal effects case and 5 vs. 5 normal effects case. The only differences are that filtering fails to beat ash (but it still beats apelgm) and more truly large effects are incorrectly and overly shrunk to zero by ash but not overly shrunk by apeglm.